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Abstract. 

A new approach is laid out to investigate the two-photon atomic transitions. It 
is based on application of the finite basis solutions constructed from the Bernstein 

Polynomial (B-Polynomial) sets. Wc show that such an approach provides a very 
promising route for the relativistic second- (and even higher-order) calculations since it 
allows for analytical evaluation of the involved matrices elements. In order to illustrate 
possible applications of the method and to verify its accuracy, detailed calculations are 
performed for the 2si/2 — > lsi/2 transition in neutral hydrogen and hydrogen-like 
ions, and are compared with the theoretical predictions based on the well-established 
B-spline-basis-set approach. 



PACS numbers: 31.15.-p, 31.30.Jv, 32.70.Fw, 32.80.Wr 



Submitted to: J. Phys. A: Math. Gen. 



20 January 2013 



Calculation of the two-photon decay rates of hydrogen-like ions 



2 



1. Introduction 

High-order perturbation calculations in atomic physics generally require summations 
over the complete spectrum of the system under consideration. Within the relativistic 
framework, such a "summation" is not a simple task since it includes a summation 
over discrete part of the Dirac spectrum as well as the integration over the positive- as 
well as negative-energy continua. A number of methods have been developed over the 
past decades to perform this summation consistently. Apart from the various Green's 
function approaches [HIS], the discrete-basis -set method is widely employed nowadays 
in (relativistic as well as non-relativistic) high-order calculations [21 IH El El C] ■ In this 
method, a finite set of discrete pseudostates is constructed from some basis functions 
and utilized for carrying out the summation. The particular choice of a suitable set 
of basis functions is crucial for the practical implication of the method. Usually, the 
discrete (pseudo-) solutions are built up from piecewise polynomial sets. The piecewise 
polynomials are precisely defined, can be calculated rapidly on modern computer 
systems, and can represent a great variety of functions. They can be differentiated 
and integrated easily [S]. 

The basis splines |9], also called B-splines, are one of the most commonly used 
family of piecewise polynomials. These polynomials, which are well adapted to numerical 
tasks, have been successfully used in many atomic-physics studies. For example, 
Johnson and co-workers have applied the B-splines to the many-body perturbation 
theory [TUl [H]. Froese Fischer and co-workers used it in (variational) Hartree-Fock 
calculations and continuum problems [121 US]- Qiu and Froese Fischer introduced 
the integration by cell algorithm for Slater integrals in a B-spline basis obtaining an 
improved efficiency and accuracy over traditional methods p3l [15]. Bhatti and co- 
workers [16] used similar techniques to find an approximate solution of a set of the non- 
homogeneous second-order differential equations and to obtain static polarizabilities 
of hydrogenic states. Indelicato and co-workers employed B-splines in the multi- 
configuration Dirac-Fock (MCDF) relativistic atomic structure calculations [T71 ITS] 
and in relativistic two-photon decay calculations [H [IHl [20] . 

While B-spline basis sets were proven to be an important tool for studying the 
variety of atomic structure and dynamics problems (see Ref. [31 |2T] for further details 
and examples), one might adopt other piecewise polynomial sets to speed up relativistic 
high-order calculations as might be highly required, for example, for studying the two- 
photon transitions in many-electron systems. In this work we argue that the Bernstein, 
or B-polynomials [2^ may serve as a good alternative to the B-splines since they 
allow for analytical finite-basis-set calculations. These are polynomial functions of n^^ 
degree that have been recently used to obtain the solution of some linear and non-linear 
differential equations [23l EH [25]. Bhatti and Perger [26] developed an algorithm for 
constructing accurate solutions to the radial Dirac equation in a B-polynomial basis 
set. By using this algorithm and the Galerkin method, accurate calculations have been 
performed for the bound-state energies of hydrogenic systems. 
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In this work, we employ the finite (discrete) solutions constructed from the B- 
polynomial sets in order to explore two-photon decay of hydrogen-like ions. Theoretical 
analysis of this process requires evaluation of the second-order transition amplitudes 
and, hence, can be used as a "testing ground" for the high-order B-polynomial 
calculations. To explain the background of these calculations we will recall in Sections 



|2.2| and |2.3| the application of finite-basis-set methods for dealing with the Dirac 
problem. In particular, we will derive the generalized eigenvalue problem whose solutions 
form a complete set of atomic pseudostates. In order to find these solutions, however, 
one has to agree first about the explicit form of the basis functions. In Section 2.4, we 



will introduce B-polynomial basis sets, and discuss their properties. We will use these 
sets to construct the (pseudo) spectrum of the ion (or atom). Summation over such a 
spectrum, which appears within the second-order perturbation theory, will be performed 
later in Section 2^ and will allow us to derive the two-photon decay rates. Results of 
the relativistic calculations for these rates obtained for the 2si/2 — ?■ lsi/2 transition in 
neutral hydrogen and hydrogen-like ions will be presented in Section |3] Apart from 
the results obtained by making use of the B-polynomial sets, we present here also the 
"standard" B-spline calculations as well as the predictions by Labzowsky et al f27\ . 
Detailed comparison with these predictions will allow us to justify the application of 
the B-polynomial aproach in second-order calculations and to underline its advantages. 
The summary of our work will be given finally in Section |4| 

The atomic system of units {e = = h = 1) is used throughout the paper unless 
otherwise stated. 



2. Theory 

2.1. General approach 

The second-order relativistic calculations in atomic and molecular physics often require 
summation over the complete spectrum of the system under consideration. Such a 
summation, that includes integration over the positive- as well as negative-energy 
continua, can be performed very efficiently if one considers the (atomic or molecular) 
system to be enclosed in a finite cavity with a radius R. This allows for a discretization 
of the continua and, hence, for a representation of the entire Dirac spectrum in terms 
of the pseudostate basis functions. A (quasi-complete) finite set of these functions are 
determined subsequently by making use of the variational Galerkin method |28j . 



2.2. Finite basis set approach to the Dirac equation 

Having briefiy discussed the general context, we are ready now to apply the finite basis 
set method for solving the eigenvalue Dirac's problem 

(ccK ■ p + (3c^ + V{r)) u{r) = eu{r) , (1) 

where a and (3 are the usual 4x4 Dirac matrices, and V{r) describes Coulomb 
interaction between an electron and nucleus. Moreover, in Eq. ([T]) we have replaced the 
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total electron energy E hj e = E — c^ to render easy the comparison with nonrelativistic 
calculations. 

Since the potential V{r) is central, the eigenfunctions of the Dirac Hamiltonian can 
be written in the standard form as 



Ur 



(r) 



riK ) 
riK 



Km \^ ) 

(r) 



(2) 



with VL^rnii) being the Dirac spin-angular function. The angular quantum number k is 
defined by 



if 
if 



-1/2 
+ 1/2 



(3) 



where I and j are the electron orbital and total angular momenta, respectively. By 
substituting the wavefunction ^ into Eq. ([T]) and by performing some simple angular 
momentum algebra one can obtain the conventional set of radial Dirac equations 



V{r) 
-c 01 



c 01 
-2meC^ + V{r) 
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to determine the large, PriK^r), and the small, Qr 
of shortness, we introduced here the operator 



O^ 



d ^ K 
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(4) 

radial components. For the sake 

(5) 



For the further evaluation of Eq. Q we assume the ion (or atom) is enclosed 
in a finite cavity with a radius R large enough to get a good approximation of the 
wavefunctions, with some suitable set of boundary conditions. In order to construct 
these functions, we shall turn to the principle of least action [10] 

SS, = 0, (6) 

from which the Dirac equation can be derived. In this expression, the action is 
defined as 

S. = l r {cPn.{r)01QUr) - cQUr)OlPn.{r) 
^ Jo 

+ Vir) [PUrf + QrUr?] -2m,c'Qn^{rf} dr 



[Pn.{rf + Qnn{rf] dr + 



bond 



(7) 



where the upper integration limit R is the radius of the confining cavity. The term S^"^"^, 
that will be specified below, stands for the boundary conditions and the parameter e is 
a Lagrange multiplier introduced to ensure the normalization constraint. 



I [Pnn{rf + Qn.{rf]dr = l. 
Jo 



(8) 
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Here, the large, P„ 



and small, Qr 



radial components of the electron 



wavefunctions can be written as a finite expansion 



N 



i=l 
N 



Q{r)= Y,<1^B^ 



(9) 



i=l 



over some basis functions Bi{r). The explicit form of these functions is not crucial for 
the following discussion and will be specified later in Section 2.4.2 In Eq. ([9]), moreover, 
the subscripts n and k have been omitted from the functions Pnnir) and QriKir) for the 
sake of notation simplicity. 

By inserting now the radial components ^ into the least action principle ([6]) and 
by evaluating the variation with respect to change of expansion coefficients pi and 
gj, we obtain the matrix equation 

Av = eBv, (10) 

to determine the vector v = {pi,p2, . . . ,Pn,(1i,(12, ■ ■ ■ ,(1n) and where A and B are 
symmetric 2N x 2N matrices given, respectively, by 

.r . 



A 



— c 



(V) 
iD) + 
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r / J 



c[iD) 
-2mec\C) + {V) 



+ A 



bond 
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(C) 
(C) 



11 



(12) 



The matrix A^°^'^ refiects the boundary conditions, and the N x N matrices (C), (D), 
{V) and (n/r) are given by 



Bi{r)Bj{r)dr, 
B^{r)^Bj{r)dr, 



B,{r)-Bj{r)dr. 



Bi{r)V{r)Bj{r)dr . 



(13) 
(14) 
(15) 
(16) 



Equation ( 10 ) is known as a generalized eigenvalue problem that can be solved by 
employing linear algebra standard techniques. In the present work, for example, we 
have used the well-established LAPACK 3.3.0 package [29j. By using this package we 
obtain 2A^ real eigenvalues eA and 2A^ orthogonal eigenvectors that span both positive 
and negative energy solutions. Solutions labeled by i = 1, describe the continuum 
e\ < —2mc^ and solutions labeled by i = A^ + 1,...,2A^ describe bound states and the 
positive continuum > 0. 
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2.3. Spurious states and boundary conditions 

The practical implementation of the finite basis set approaches is usually complicated 
by the well-known problem of spurious states. These non-physical states appear as 
solutions of the single-particle radial Dirac equation for k > (pi/2, ^3/2, ■ ■ ■ orbitals) 
[30] . Although the spurious solutions "spoil" the spectrum of the ion (or atom) under 
consideration, they are required for providing completeness of the basis set. The problem 
of spurious states has been discussed in detail in the literature, and several solutions 
were proposed [301 EB EH |33l [31] • Johnson an co-workers [10], in their pioneering 
applications of the B-splines to the relativistic many-body problem, have suggested to 
adopt the function S^°^^ in Eq. ^ as 

' i[p\R)-Q\R)]-ipmp{^)-Qm 

cBond _ J for /t < , , 

^ l[P\R)-Q\R)]-lP{Q)[cP{0)-Q{Q)] ' ^ > 

for K > 

in order to lift the spurious states to lower energies (to the negative continumm), thus 
restoring the low-energy mapping to the physical solutions. For variations of P(r) and 
Q{r) the boundary terms vanish if 

P(0) = and P(P) = Q{R). (18) 

The latter boundary condition at the outer boundary r = P is the MIT-bag-model 
condition [35], and was included to avoid problems with Klein's paradox that arises 
when one attempts to confine a particle into a cavity, essentially by forcing the radial 
current crossing the boundary to vanish 

2.4. Finite- basis expansion 



2.4.1. B-polynomials expansion As seen from Eqs. (10)-(16), any further analysis of 
the generalized eigenvalue problem requires the knowledge on the the explicit form of 
the basis functions Pj(r). In the present work we will construct these functions from 
the B-polynomial as well as B-spline spline sets. While the latter case will be discussed 



in Subsection |2.4.2[ here we shall recall the main features of B-polynomials, also known 
as Bernstein functions. 

The B-polynomials of fcth-order are defined by [221 ESI El] 

B.Ar)={ .-0.1.....^. (19) 

where the standard form of the binomial coefficients are utilized, 
k\ k\ 

and where a and h denote the limits of the interval [a, h] over which the polynomials are 
defined to form a complete basis. Since the atomic system is defined in a finite cavity 
of radius P, we take a = and b = R. 
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As seen from definition (19), there are are {k + 1) polynomials of degree k. By 
definition, we set Bi^k{r) = if z < and i > k. As an example, a set of 11 B- 
polynomials of degree 10 is plotted in Fig. [l] where where it is shown that each B- 
polynomial is positive and the sum of all B-polynomials is unity. 

The great advantage of B-polynomials is that they allow for an analytical evaluation 
of the matrices (C), (D) and (ft/r) in the generalized eigenvalue problem (10). That is, 
by inserting Bi = Bi^k{r) into Eqs. (13)-(15), we obtain 

k 



r/ij 




(21) 



(22) 



(23) 



Apart from the basis functions Bi^k{r), the knowledge on the nuclear charge distribution 
is also required for an evaluation of the matrix (V") whose elements are defined by 
Eq. (16). For the point-like nucleus potential, V^{r) = —Z/r, for example, these 
matrix elements read as 




(24) 



A more complicated expression for the matrix {V) is obtained to account for the finite 
nuclear size effects. To address these effects, we employ here the potential 

2 

r < tn 



V'^(r) 



z 

2rN 



(25) 



-7, r>rN 
due to a uniform spherical nuclear charge distribution with radius r^. By employing 
this potential in Eq. (16), we finally obtain 



(\/p 



r^Y+^ 82^1 [{l + t+J,t+J -2k},{2 + t + j}/^] 



+ - 



R 
R 



1 + 1+ J 

+j 2F, [{3 + z+j,t+j -2k},{A + i + j}/-^] 



(26) 



3 + t + j 

where B {x; h, k) is the incomplete beta function, pFq is the generalized hypergeometric 
function, and ()^ is the Pochhammer symbol. 
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2.4-2. B-splines expansion Since the B-spline basis set approach has been discussed 
in detail elsewhere [8l |10] , here we will restrict ourselves to a very brief compilation of 
its basic expressions. Following de Boor [9], we divide the interval of interest [0, R] into 
segments whose endpoints define a knot sequence {ti} = 1,2, . . . ,n + k. The B-splines 
of the order k, Bi^k{r), are defined on this knot sequence by the recurrence relation 

B^,k{r) = -^—^B,^k-i{r) + ~/ ^m.fc-i(r) , (27) 
where the B-splines of the first order read as 

Note, that in these expressions the number of knots ti is by k larger than the 
number of splines. The first and the last k knots must be equal and are defined as: 
^1 = ^2 = • • • = ^fc = and = t„+2 = . . . = tn+k = R- 



2.5. Two-photon decay rates 

In the previous Sections we have obtained the finite (discrete) basis solutions of the Dirac 
eigenproblem ([T]) constructed from the B-polynomial as well as B-spline sets. Now we 
are ready to apply these solutions for studying the two-photon transitions in hydrogen- 
like ions. Not much has to be said about the theoretical background for describing 
such a second-order process. In the past, relativistic calculations of both, the total 
and the differential two-photon decay rates have been discussed in detail elsewhere 
[U Uni EDI EZ]- Below, therefore, we will repeat just basic expressions, relevant for 
discussing the role of the finite basis-sets in the two-photon calculations. 

Usually, the properties of the two-photon atomic transitions are evaluated within 
the framework of the second-order perturbation theory. This theory provides the 
following expression for the differential in energy decay rate |38j . 

-{f\Al\v){v\A\\i) , {f\A\\v){v\Al\i)^ ' 



dw UJ1UJ2 
duj~i ~ (27r)3c2 



E/ \J l^2l V , 

I 7? _ p. _L ,,,, 



dnidn2 , (29) 



Ey — Ei + oji Ey — Ei + U2 

where Uj is the frequency and dQj is the element of solid angle of the j^^ photon. In this 
expression, moreover, \i) = \niji^i), |z/) = \{eu)nuji,fiy) \f) = \nfjffif) denote solutions 
of the Diracs equation for the initial, intermediate and final ionic states respectively, 
while Ei, Ey and Ef are the corresponding one-particle energies. Because of energy 
conservation, Ei and Ej are related to the energies uji^ of the emitted photons by 

Ef - Ei = uji+ uj2 = u}^. (30) 

For a photon plane-wave with propagation vector kj and polarization vector e.j 



{e.j ■ kj=0), the operators A* in Eq. (29) are given by 



A* = a- (e, + Gk^ e-'^^r _ Ge'^'^rr (31) 
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where ck is the vector of Dirac matrices and G is an arbitrary gauge parameter. Among 
the large variety of possible gauges, Grant [39j showed that there are two values 
of G which are of particular utility because they lead to well-known nonrelativistic 
operators. If G = 0, one has the so-called Coulomb gauge, or velocity gauge, which 

1 /2 

leads to the dipole velocity form in the nonrelativistic limit. If G = [(L + 1) /L] ' , for 
example, G = a/2 for El transitions {L = 1), one obtains a nonrelativistic expression 
which reduces to the dipole length form of the transition operator. From the general 
requirement of gauge invariance the final results must be independent of G. The gauge 
invariance of the two-photon relativistic calculations was studied by Goldman and Drake 
[S7] and by Santos et al 

So far, we have discussed the general expressions for the two-photon transition 
rates that are differential in the energy ui of one of the photons. By performing an 
integration over this energy one may easily obtain the total rate 

/dw 
dco~ 



which is directly related to the lifetime of a particular excited state against the two- 
photon decay. In this expression, we introduce the transition energy as Ut = Ui + U2- 

The evaluation of the matrix elements (/ \A*\ z/) and (z/ \A*\ and, hence, of the 
differential (29) as well as total (32) rates involves the radial integrals 



and 



u) = [Pf{r)Q,{r) ± Qf{r)P,{r)] jl (^) rfr, (33) 

Jl{uj) = [Pfimir) + Q/(r)g.(r)] (^) dr, (34) 

where ji (x) is the spherical Bessel function of the first kind ^1] . By making use of a a 
finite piecewise polynomials basis-set to describe the large and small radial components 
of the Dirac wavefunctions [cf. Eq. (P)], both, /^(w) and Jl(ci;) can be reduced to a 



linear combination of the integrals 

F'^'ir)F^^Hr)jL (^) dr = J2Yl ft'^f^'^ij^h ■ (35) 

i=i j=i 

Here, for the sake of shortness, we denote F^{r) = P{r) and F^{r) = Q{r) as well as 
= Pi and /f = qi, and the matrix elements {jL)ij are defined by 

= B,{r)B,{r)j, (^) dr . (36) 

These radial matrix elements are the "building blocks" used to evaluate two-photon 
decay rates. An efficient and fast evaluation of these matrix elements is crucial for 
studying transitions not only in hydrogen-like but also in many-electron ions. Such 
an evaluation can be easily performed if we employ the B-polynomials as basis-set 
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functions. Inserting Eq. (19) into (36) we find that {jL)ij matrix elements are given 
analytically by 

^ ' k \ ( k \ TT{i+j + L)\{2k-i-j)\ 

i 



L)ij 



R 



uR 



J I 22(-^+^+i) 
i+j + L + 1 i+j + L + 2' 



2 ' 2 
2L + 3 2A; + L + 2 2k + L + 3 



uR 
~2^ 



where pF„ is the regularized generalized hypergeometric function, 



,Fg [{ai, ...,ap} , ,x] 



1 



and r(x) is the gamma function. 



3. Results and discussion 



(37) 



E^H^T^^' (38) 



3.1. Determination of the optimal set of parameters 

Having discussed the application of the finite basis solutions constructed from the B- 
spline and B-polynomials sets for the two-photon relativistic studies, we will employ 
now these solutions to analyze the properties of the 2s 1/2 — >■ lsi/2 transition in the 
hydrogen-like ions. Before we start an analysis, however, we must find the optimal set 
of parameters^ such as the number of the basis functions or the size R of the cavity 
radius, to be used in these calculations. 

To determine the optimal set for the B-polynomial approach, we consider first the 
transition energy [Eq. (30)]. Accurate computation of this energy is crucial since 



both the differential [Eq. (29)] and total [Eq. (32)] decay rates depend quadratically 
on it. Detailed calculations of this energy value were carried out by making use of 
B-polynomial basis sets for various numbers of polynomials, ^bp, and of the cavity 
radius, R. For the particular case of the point-like nucleus, results of our B-polynomial 
calculations have been compared against the predictions of the well-known expression 
for the Dirac bound-state energies (see, e.g., Ref. ^2] ) 

-I -1/2 



_ Exact 



1+ T 



n 



■n 2 



(39) 



where n is the principal quantum number and n is the angular quantum momentum of 
electron [cf. Eq. (l3l) 



Apart from the energy values , evaluation of two-photon transition rates (29) 



and (32) requires also detailed knowledge of the atomic wavefunctions. As usual in 



atomic structure calculations, an indication for the completeness and quality of the basis 
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set can be obtained from the comparison of the results obtained within two different 
gauges [cf. Eq. (31)]. In this way, detailed calculations of the total decay rates for the 
leading two-photon 2E1 2si/2 — ?■ lsi/2 channel have been performed, in both length 
and velocity gauges. The gauge invariance condition together with the the energy ojt 
accuracy determine the optimal set of parameters. 

The relative difference A^^ = Kuf"^^^^ — ujf^)/oof'^^'^^\ between the exact cjj^^'*^* 
solution and the basis set value uf^ is presented in Fig. |2] as a function of tt-bp for 
Z equal to 1, 40 and 92. We display in Fig. |3]the relative difference between the length 
and velocity gauge decay rate values of the 2E1 2si/2 —J- lsi/2, '^i-v, as function of the 
number of B- Polynomials, nep, for the same values of the atomic number Z. 

As seen from Figs. |2] and [3} the optimal wbp value for Z = 1 is greater than 34 
and the optimal R is 50. For Z = 40 the parameters are tt-bp > 40 and R = 1, and for 
Z = 92 are hbp > 40 and R = 0.25. 

Working in double precision, we noticed a loss of numerical significance in the 
results for tt-bp > 23 due to the LAPACK 3.3.0 routines used for the solution of the 
generalized eigenvalue problem. We detected that this is due to the difficulty of these 
routines, namely the DSYGVX routine, to deal with diagonally dominant matrices that 
have very large diagonal elements and very small off diagonal elements. This problem 
was overcome by using quadruple precision for the evaluation of the matrix elements and 
compiling the LAPACK subroutine in quadruple precision as well, using the compiler 
auto-doubling option. This situation is illustrated in Fig. |4| where it is plotted the 
relative difference Ai_v for Z = 1 and R = 40, obtained in double and quadruple 
precision. 

Nevertheless, we should emphasize that even with a low nBp value, such as 20, and 
working in double precision, we get a gauge invariance in the two-photon rate less than 
10-'^. 

Before turning to the second-order relativistic calculations performed within the 
finite-basis-set approaches, let us recall that apart from the B-polynomials, we used 
also the "standard" B-spline solutions. For these solutions, we adopted the following 
parameters described in Ref. [19]: k = 9, n-ss = 60, and R = 60 a.u. Moreover, the 



integration over the photon frequency in Eq. (32) has been performed using a 15-point 
Gauss-Legendre algorithm. 



3.2. Second-order relativistic calculations 

The "optimal" sets of parameters for both, the B-spline and B-polynomial approaches, 
determined in the previous section will be used below for the computation of the two- 
photon transition rates in quadruple precision. However, first let us briefly return to the 
energy calculations performed in these two approaches. In Table [T]we report the relative 
differences between the computed energies of nsi/2 states of neutral hydrogen and the 
corresponding exact values e^^^'^^ defined by Eq. 
represent the differences obtained by employing the B-polynomial and B-spline sets. 



39). Here, AE^rf°' and AE^-^p^'"*^^ 
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correspondingly. Moreover, AE^p^°^ denotes the relative difference calculated with the 
B-polynomial values by Bhatti and Perger [21] . The B-polynomial relative differences 
are very similar with the B-Spline relative differences, for n > 2. For n = 1,2, the 
former method achieved an excellent 12 digits agreement with the exact result, two 
orders of magnitude better than the B-Splines. 

In Table [2] we display the most significant multipole contributions to the 2si/2 — ^ 



lsi/2 two-photon total decay rate [Eq. (32)] of neutral hydrogen. Again, calculations 
have been performed within two different gauges and by employing B-polynomial as 
well as B-spline basis sets. Multipoles with only magnetic components such as 2M1 
only have velocity gauge. As seen in this Table, both approaches yield almost identical 
results (with the relative error of less than 2 x 10~^) for all multipole decay channels. 
Furthermore, one may observe a remarkable agreement (smaller than 10~^^) between 
the values obtained in the length and velocity gauges for both basis sets. 

It should be emphasized that with the B-Polynomials approach, that employs 
analytical evaluation of the second-order matrix elements [cf. Eqs. (33)- (37)] and 
uses a smaller optimal nep, the computation time required for the determination of 
each multipole contribution is, in quadruple precision, about two and half times smaller 
than the time required by the B-Splines methods. These results clearly indicate that 
the finite-basis-set approach based on B-polynomial solutions provides an alternative 
tool for studying the two-photon transitions. 

Besides the neutral hydrogen, extensive test of the B-polynomial approach has been 
performed also for the two-photon transitions in medium- and high-Z hydrogen-like 
ions. In Table[3]we display, for example, the total rates of the leading, 2£^1 2si/2 — )• lsi/2 
channel for selected values of the nuclear charge Z. Here we used both, B-spline and B- 
polynomial approaches to carry out the intermediate-state summation in Eq. (29) over 
the complete Dirac's spectrum {W^) as well as over the positive- (W^) and negative- 
energy {W~) solutions only. The role of Dirac's continua in relativistic second-order 
calculations has been the subject of recent theoretical investigations [201 EZ] partially 
because of their impact on (future) many-body studies. As seen from Table [3] , the 
good agreement between the predictions obtained within the B-polynomial and B-spline 
approaches can be found for the ions along the entire isoelectronic sequence. Moreover, 
both our calculations show a perfect agreement between gauges, which is better than 
the one reported by Labzowsky et al [27] . 

Again, the presented calculations validate the B-Polynomial basis set as a useful 
and appropriate basis set to study two-photon decays in atomic systems, and various 
atomic calculations that involve a summation over the Dirac spectrum. 



4. Conclusions 



In the present work we have investigated, for the first time, the efficiency and accuracy 
of the B-Polynomial basis set for studying the two-photon transitions in hydrogen- 
like systems. Based on the finite-basis-set approach, the generalized eigenproblem in 
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Eq. ( 10 ) was solved in order to provide the eigenvalues and eigenfunctions, which were 
successfully used for the calculations of the two-photon decay rates. We took advantage 
of the B-Polynomial properties and derived fully relativistic analytical expressions for 
the two-photon rates both, within the pointlike- and finite-nucleus models, instead of 
employing numerical methods. 

In order to illustrate the application of the B-polynomial method and to verify its 
accuracy, we have performed detailed calculations of the total rates for the 2si/2 
two-photon transition in neutral hydrogen and hydrogen-like ions. Results of these 
calculations have been compared with the predictions of the well-established B-spline 
approach. While the perfect agreement between the results of both, B-spline and B- 
polynomial approximations was observed along the entire isoelectronic sequence, the 
B-polynomial calculations were found to be much less computationally demanding. 

It was noticed that if we consider a basis set with more than 23 B-Polynomials, 
we need to work in quadruple precision due to the LAPACK 3.3.0 routines limitation 
in dealing with matrices that have very large diagonal matrix elements and very small 
off-diagonal values. 

We conclude that the B-Polynomial basis sets are suitable for investigation of the 
two-photon transitions with the great advantage of enabling analytical expressions for 
the involved integrals, which speed up the calculations. 

The presented procedure may be easily extended to calculate two-photon 
transitions in many-electron ions or atoms, as well as to the quantum electrodynamics 
(QED) calculations that usually require the summation over the complete Dirac 
spectrum. 
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Table 1. Relative differences between the computed energy eigenvalues of hydrogen 
(s states) using B-polynomials and B-splines and the exact results of the Coulomb- 
field Dirac equation, E^''^'=\ given by Eq. AE^^^°^ and Ai;^^ represent, 
respectively, the relative differences calculated in this work with the B-polynomials 
and with the B-splines. AE'gp denote the relative difference obtained with the B- 
polynomials values by Bhatti and Perger [24] . Powers of ten are given in parentheses. 

Exact A pB-Pol a T-iB-splincs a pB-Pol 



1 


-0.5000066565953603 


1.6(- 


-12) 


3.8(- 


-10) 


2.7(- 


-12) 


2 


-0.12500208018900594 


1.4(- 


-12) 


1.4(- 


-10) 


1.5(- 


-8) 


3 


-0.05555629517766647 


1.7(- 


-7) 


3.4(- 


-9) 


1.8(- 


^3) 


4 


-0.03125033803007682 


1.4(- 


-3) 


5.9(- 


'5) 


2.3(- 


-1) 



Table 2. Multipole contributions (in s""'^) of the 2si/2 ^ 1*1/2 two-photon decay 
for Z = I. The relativistic calculations have been performed within the velocity and 
length gauges, using the B-polynomials and B-splines basis sets. Ai_v stands for the 
relative difference between the length and velocity gauge values. Powers of ten are 
given in parentheses. 

Multipoles Contribution (s^^) 



B-Polynomials B-Splines 



lenght gauge Ai_v length gauge Ai_v 



2£;i 


8.2290591586 




< 1.0(- 


-26) 


8.2290591509 




< 1.0(- 


-15) 


E1M2 


2.5371807735(- 


-10) 


< 1.0(- 


-25) 


2.5371807635(- 


-10) 


< 1.0(- 


-15) 


2M1 


1.3803580496(- 


-11) 






1.3803580473(- 


-11) 






2E2 


4.9072289232(- 


-12) 


< 1.0(- 


-34) 


4.9072289165(- 


-12) 


< 1.0(- 


-14) 


2M2 


3.0693510074(- 


-22) 






3.0693509833(- 


-22) 






E2M1 


1.6393565197(- 


-23) 


< 1.0(- 


-34) 


1.6397413530(- 


-23) 


< 1.0(- 


-4) 


Total 


8.2290591589 








8.2290591512 
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Table 3. Total two-photon 2E1 2si/2 — > lsi/2 decay rates (in s^^) for selected values 
of the nuclear charge Z. The relativistic calculations have been performed within the 
velocity and length gauges, using the B-polynomials and B-splines basis sets, and by 
carrying out intermediate-state summation over the complete Dirac's spectrum (M^"^) 
as well as over the positive- (W^) and negative-energy (W^) solutions only. Ai_v 
stands for the relative difference between the length and velocity gauge values, and 
LSS denotes the values calculated by Labzowsky et al [57]. Powers of ten are given 
in parentheses. 



Z Decay rates (s ^) 



B-Polynomials B-Splines LSS 







Icnght gauge 




■V 


length gauge 


Ai- 


V 


length gauge 


Ai_v 




1 


W+ 


8.22861 






8.22861 






8.2206 








w- 


2.49477(-8) 






2.49481(-8) 






3.9975(-22) 










8.22906 


< 1.0(- 


-26) 


8.22906 


< 1.0(- 


-15) 


8.2207 


< 1.0(- 


-4) 


40 


W+ 


2.96130(-l-10) 






2.96124(-hl0) 






3.1953(-t-10) 








w- 


2.10270(-h8) 






2.10250(+8) 






5.8284 










3.19862(-fl0) 


< 1.0(^ 


-13) 


3.19858(+10) 


< 1.0(- 


-15) 


3.1953(-M0) 


< 1.0(- 


-4) 


92 


W+ 


2.90482(-|-12) 






2.90409(+12) 






3.8230(-hl2) 








w- 


6.85553(-l-ll) 






6.80648(-t-ll) 






1.2851(-H5) 










3.825839(-l-12) 


< 1.0(- 


-9) 


3.82555(+12) 


< 1.0( 


-15) 


3.8216(4-12) 


< 1.0(- 


-5) 
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Figure 1. The set of 11 B-Polynomials of degree 10 are shown in the interval 
The quantities are dimensionless. 
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Figure 2. Relative difference A^^t between the computed Wt energy, given by Eq. [30) 



and the value 



.Exact 



= E. 



Exact 



Exact 



as function of the number of B-Polynomials, 



TiBp, for Z equal to 1, 40 and 92. Here, iJ^xact g^j.^ ^j^g g^act solutions of the Coulomb- 
field Dirac equation given by Eq. (39). 




Figure 3. Relative difference between the length and velocity gauge decay-rate values 
of the 2E1 2si/2 ^ lsi/2: Ai_v, as function of the number of B-Polynomials, rtep, for 
Z equal to 1, 40 and 92. 
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Z=1 and R=40 
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Figure 4. Relative difference between the length and velocity gauge decay-rate of the 
2E1 2si/2 lsi/2, Ai_v, for Z = 1 and R = 40, obtained in double and quadruple 
precision. 



